function ana_out = lmj_modeanalysis_compute(funcmod, pertmat, data, in, varargin)
% XIPLOT_ALL
% Given multiple group of parameter values
% evaluate the model and plot autocovriance, variance decomposition, IRF,
% impulse responses.
% apply Kalman Filter to data and plot states, errors and
% state-decomposition
%
% INPUT:
%       "FUNCMOD":      function handle of the model
%       "PERTMAT":      [NPar NPert] matrix of parameters.
%       "DATA":         [NOBS NSER] matrix of data
%       "IN":           Structure that contains other necessary input
%          .STATE_SEL:  Cell array of names of selected state
%          .STANAMES:   names of all STATES
%          .SHONAMES:   names of all shocks
%          .SAVEPATH:   name of the save path
%          .PARNAMES:   name of parameters
%          .DATEVEC:    vector of dates
%          .MODELNAMES: names of the models
%          .FLAG_BLOCK: Whether to put
% OUTPUT:  LOTS and LOTS of graphs
%
%
% Dependencies (mfiles and matlab functions)
%   (major functions)
%       kfilter_all.m
%       momentszselstate.m
%       plot_multiple3d.m
%       accov.m
%       priorposterior_acout.m
%       plot_decompositions.m
%       multimod_acout.m (subfunction)
%
%   (tools) strdate, cr_dir, cellposition, dispaj, suptitle, ps2pdf, vline
%
% For: Alejandro Justiniano
% Research Assistant: Xi Luo
% Created: 02/22/2010

cucd = pwd;

% Parse INPUTS
try 
    state_sel = in.state_sel; % commented out 07/06/2010
catch
    disp('no particular states are selected')
end 
stanames =  in.stanames;
shonames =  in.shonames;
save_path = in.savepath; 
modelnames = in.modelnames;


[numpar, numpert] = size(pertmat);
try
%[state_pos,ana_in.flag_adds]=extractlevels(state_sel,stanames); % commented out 07/06/2010
[state_pos,in.flag_adds]=extractlevels(state_sel,stanames); %MO 10/24/2011
catch
    display('Error extracting state positions');  %MO 10/24/2011 
    error('Error extracting state positions'); 
end


%state_pos = cellposition(state_sel, stanames);

% Solve for each PARAMETER.
numst    = length(stanames);
numshock = length(shonames);
pertv_ok = zeros(numpert, 1);
GGmat    = zeros(numst, numst, numpert);
CONSmat  = zeros(numst, 1, numpert);
RRmat = zeros(numst, numshock, numpert);
SDXmat = zeros(numshock, numshock, numpert);
eumat = zeros(2, 1, numpert);
%ZZmat = zeros(size(data,2),numst,numpert ); 

for ii = 1:numpert
    [GGtemp,RRtemp,CONStemp,eutemp,SDXtemp,ZZtemp, outtemp]=feval(funcmod,pertmat(:, ii), varargin{:});
    if isequal(eutemp,[1;1])
        dispaj('pertvec Number ',ii,'  will be plotted');
        pertv_ok(ii)=1;
        GGmat(:, :, ii) = GGtemp;
        RRmat(:, :, ii) = RRtemp;
        CONSmat(:, :, ii) = CONStemp;
        eumat(:, :, ii) = eutemp;
        SDXmat(:, :, ii) = SDXtemp;
        ZZmat(:, :, ii) = ZZtemp;
        
        try
            AA = outtemp.AA; 
        catch 
            AA = []; 
        end
        if ~isempty(AA); 
            AAmat(:, :, ii) = AA; 
        end
        
    else
        dispaj('pertvec Number ',ii,'  cannot be plotted');
        pertv_ok(ii)=0;
    end
end
clear GGtemp RRtemp CONStemp eutemp SDXtemp ZZtemp ii

posok = find(pertv_ok == 1);
if length(posok(:)) < 1
    error('No parameter in PERTMAT delivers unique and stable solutions')
end
disp([num2str(length(posok)), ' out of ', num2str(numpert), 'yields valid results']);

pertmat = pertmat(:, posok); % matrices of parameters that yield solution
GGmat   = GGmat(:, :, posok);
RRmat   = RRmat(:, :, posok);
CONSmat = CONSmat(:, :, posok);
eumat   = eumat(:, :, posok);
SDXmat  = SDXmat(:, :, posok);
ZZmat   = ZZmat(:, :, posok);

if exist('AAmat', 'var');
    AAmat = AAmat(:, :, posok);
end

numpert = size(pertmat, 2);
modelnames = modelnames(posok);

if ~exist('state_pos', 'var')
    state_pos = []; 
end 
%% Starting the BIG LOOP HERE

for ii = 1:numpert
    dispaj('begin parameter', ii);
    GG = GGmat(:, :, ii);
    RR = RRmat(:, :, ii);
    CONS = CONSmat(:, :, ii);
    SDX = SDXmat(:, :, ii);
    ZZ = ZZmat(:, :, ii);
    if in.flag_mixfreq == 1
        AA = AAmat(:, :, ii); 
        %out = momentszselstate_mixfreq(GG,CONS,RR,SDX,ZZ,AA, in.case_tag,in,state_pos);
        try 
            numit = outtemp.NUMIT; 
        catch 
            % assuming monthly to quaterly aggregation
            numit = 3; 
        end
        out = momentszselstate_mixfreq_general(GG,CONS,RR,SDX,ZZ,AA, in.case_tag,in,state_pos, numit);
    else 
        out = momentszselstate(GG,CONS,RR,SDX,ZZ,in,state_pos);
    end
        
    % STORing in MATRICES
    stdthz_mat(:, :, ii) = out.stdthz;
    stdths_mat(:, :, ii) = out.stdths;
    acthz_mat(:, :, :, ii) = out.acthz;
    acths_mat(:, :, :, ii) = out.acths;
    covthz_mat(:, :, :, ii) = out.covthz;
    covths_mat(:, :, :, ii) = out.covths;
    vdecstaz_mat(:, :, ii) = out.vdecstaz;
    vdecstas_mat(:, :, ii) = out.vdecstas;
    irfz_mat(:, :, :, ii) = out.irfz;
    irfs_mat(:, :, :, ii) = out.irfs;
    vdechz_mat(:, :, :, ii) = out.vdechz;
    %vdechs_mat(:, :, :, ii) = out.vdechs; %commented out 07/06/2010
       
    stdsimz_mat(:, :, ii) = out.stdsimz;
    if ~isempty(out.stdsims)  % if non empty
    stdsims_mat(:, :, ii) = out.stdsims;
    end 
    acsimz_mat(:, :, :, :, ii) = out.acsimz;
    acsims_mat(:, :, :, :, ii) = out.acsims;
    covsimz_mat(:, :, :, :, ii) = out.covsimz;
    covsims_mat(:, :, :, :, ii) = out.covsims;
    
    if in.flag_mixfreq
        stdsimz_tag_mat(:, :, ii) = out.stdsimz_tag;
        acsimz_tag_mat(:, :, :, :, ii) = out.acsimz_tag; 
        acsimz_mix_mat(:, :, :, :, ii) = out.acsimz_mix; 
    end 
    
    clear GG RR CONS SDX ZZ
end

    
ana_out.stdthz_mat = stdthz_mat; 
ana_out.stdths_mat = stdths_mat; 
ana_out.acthz_mat  = acthz_mat; 
ana_out.acths_mat  = acths_mat;
ana_out.covthz_mat = covthz_mat;
ana_out.covths_mat = covths_mat;
ana_out.vdecstaz_mat = vdecstaz_mat; 
ana_out.vdecstas_mat = vdecstas_mat; 
ana_out.irfz_mat = irfz_mat; 
ana_out.irfs_mat = irfs_mat; 
ana_out.vdechz_mat = vdechz_mat; 
ana_out.vdechs_mat = []; % set to be empty matrix for now.  
ana_out.stdsimz_mat = stdsimz_mat; 
ana_out.stdsims_mat = stdsims_mat; 
ana_out.acsimz_mat = acsimz_mat; 
ana_out.acsims_mat = acsims_mat;
ana_out.covsimz_mat = covsimz_mat; 
ana_out.covsims_mat = covsims_mat; 
ana_out.num_valid = length(posok); 
ana_out.pos_valid = posok; 
ana_out.modelnames = modelnames; 
ana_out.pertmat = pertmat; 

if in.flag_mixfreq
    ana_out.stdsimz_tag_mat =  stdsimz_tag_mat;
    ana_out.acsimz_tag_mat  =  acsimz_tag_mat;
    ana_out.acsimz_mix_mat  =  acsimz_mix_mat; 
end

end